Go back to the Contents page.
Press Show to reveal the code chunks.
# Set seed for reproducibility
set.seed(1982)
# Set global options for all code chunks
knitr::opts_chunk$set(
# Disable messages printed by R code chunks
message = FALSE,
# Disable warnings printed by R code chunks
warning = FALSE,
# Show R code within code chunks in output
echo = TRUE,
# Include both R code and its results in output
include = TRUE,
# Evaluate R code chunks
eval = TRUE,
# Enable caching of R code chunks for faster rendering
cache = FALSE,
# Align figures in the center of the output
fig.align = "center",
# Enable retina display for high-resolution figures
retina = 2,
# Show errors in the output instead of stopping rendering
error = TRUE,
# Do not collapse code and output into a single block
collapse = FALSE
)
# Start the figure counter
fig_count <- 0
# Define the captioner function
captioner <- function(caption) {
fig_count <<- fig_count + 1
paste0("Figure ", fig_count, ": ", caption)
}capture.output(
knitr::purl(here::here("functionality.Rmd"), output = here::here("functionality.R")),
file = here::here("purl_log.txt")
)
source(here::here("functionality.R"))Let each edge \(e\in\mathcal{E}\) be subdivided into \(n_{e}\geq 2\) regular segments of length \(h_{e}\), and be delimited by the nodes \(0 = x_0^{e},x_1^{e},\dots,x_{n_{e}-1}^{e}, x_{n_{e}}^{e} = \ell_{e}\). For each \(j = 1,\dots,n_{e}-1\), we consider the following standard hat basis functions \[\begin{equation*} \varphi_j^{e}(x)=\begin{cases} 1-\dfrac{|x_j^{e}-x|}{h_{e}},&\text{ if }x_{j-1}^{e}\leq x\leq x_{j+1}^{e},\\ 0,&\text{ otherwise}. \end{cases} \end{equation*}\] For each \(e\in\mathcal{E}\), the set of hat functions \(\left\{\varphi_1^{e},\dots,\varphi_{n_{e}-1}^{e}\right\}\) is a basis for the space \[\begin{equation*} V_{h_{e}} = \left\{w\in H_0^1(e)\;\Big|\;\forall j = 0,1,\dots,n_{e}-1:w|_{[x_j^{e}, x_{j+1}^{e}]}\in\mathbb{P}^1\right\}, \end{equation*}\] where \(\mathbb{P}^1\) is the space of linear functions on \([0,\ell_{e}]\). For each vertex \(v\in\mathcal{V}\), we define \[\begin{equation*} \mathcal{N}_v = \left\{\bigcup_{e\in\left\{e\in\mathcal{E}_v: v = x_0^e\right\}}[v,x_1^e]\right\}\bigcup\left\{\bigcup_{e\in\left\{e\in\mathcal{E}_v: v = x^e_{n_e}\right\}}[x^e_{n_e-1},v]\right\}, \end{equation*}\] which is a star-shaped set with center at \(v\) and rays made of the segments contiguous to \(v\). On \(\mathcal{N}_v\), we define the hat functions as \[\begin{equation*} \phi_v(x)=\begin{cases} 1-\dfrac{|x_v^{e}-x|}{h_{e}},&\text{ if }x\in\mathcal{N}_v\cap e \text{ and }e\in\mathcal{E}_v,\\ 0,&\text{ otherwise}, \end{cases} \end{equation*}\] where \(x_v^e\) is either \(x_0^e\) or \(x_{n_e}^e\) depending on the edge direction and its parameterization. See Arioli and Benzi (2018) for more details. Figure 1 shows an illustration of the basis function system \(\{\varphi_j^e, \phi_v\}\) on the tadpole graph (in black). Standard hat functions associated with internal edge nodes are shown in blue, while special vertex-centered function are highlighted in red.
graph <- gets.graph.tadpole(h = 1/4)
graph_to_get_loc <- gets.graph.tadpole(h = 1/40)
loc <- graph_to_get_loc$get_mesh_locations()
A <- graph$fem_basis(loc)
A <- cbind(A, A[, 1]*0) # this plots the graph
x <- graph_to_get_loc$mesh$V[, 1]
y <- graph_to_get_loc$mesh$V[, 2]
x <- plotting.order(x, graph_to_get_loc)
y <- plotting.order(y, graph_to_get_loc)
x_range <- range(x)
y_range <- range(y)
z_range <- c(0,1)
# Start plot
p <- plot_ly()
# Add function lines
for (i in 1:ncol(A)) {
z <- plotting.order(A[, i], graph_to_get_loc)
mycolor <- if (i == ncol(A)) "black" else "blue"
if (i %in% c(1,2)) mycolor <- "red"
# Add basis function trace
p <- add_trace(p, x = x, y = y, z = z, type = "scatter3d", mode = "lines",
line = list(color = mycolor, width = 3), showlegend = FALSE)
# Add vertical segments from (x, y, 0) to (x, y, z)
for (j in 1:length(x)) {
p <- add_trace(p,
x = rep(x[j], 2),
y = rep(y[j], 2),
z = c(0, z[j]),
type = "scatter3d",
mode = "lines",
line = list(color = "gray", width = 1),
showlegend = FALSE
)
}
}
z <- plotting.order(A[, 1], graph_to_get_loc)
z_new <- rep(NA, length(z))
z_new[1:11] <- 0
p <- add_trace(p, x = x, y = y, z = z_new, type = "scatter3d", mode = "lines",
line = list(color = "green", width = 3), showlegend = FALSE)
z <- plotting.order(A[, 2], graph_to_get_loc)
z_new <- rep(NA, length(z))
z_new[c(31:51, 111:121)] <- 0
p <- add_trace(p, x = x, y = y, z = z_new, type = "scatter3d", mode = "lines",
line = list(color = "green", width = 3), showlegend = FALSE)
# Display plot
p %>% layout(
scene = global.scene.setter(x_range, y_range, z_range, z_aspectratio = 1))Figure 1: Illustration of the system of basis functions \(\{\varphi_j^e, \phi_v\}\) on the tadpole graph. Notice that the sets \(\mathcal{N}_{v}\) are depicted in green and their corresponding basis functions are shown in red.
Having introduced the system of basis functions \(\{\varphi_j^e, \phi_v\}\), which will be referred to as \(\{\psi_j\}_{j=1}^{N_h}\), we can now define the finite element space \(V_h\subset H^1(\Gamma)\) as \(V_h = \left(\bigoplus_{e\in\mathcal{E}} V_{h_e}\right)\bigoplus V_v\), where \(V_v = \text{span}\left(\{\phi_v:v\in\mathcal{V}\}\right)\) and \(\dim\left(V_h\right)\) is given by \(N_h = |\mathcal{V}| + \sum_{e\in\mathcal{E}}n_e\).
Let \(\Gamma_T = (\mathcal{V},\mathcal{E})\) characterize the tadpole graph with \(\mathcal{V}= \{v_1,v_2\}\) and \(\mathcal{E}= \{e_1,e_2\}\). The left edge \(e_1\) has length 1 and the circular edge \(e_2\) has length 2. As discussed before, a point on \(e_1\) is parameterized via \(s=\left(e_1, t\right)\) for \(t \in[0,1]\) and a point on \(e_2\) via \(s=\left(e_2, t\right)\) for \(t\in[0,2]\). One can verify that \(-\Delta_\Gamma\) has eigenvalues \(0,\left\{(i \pi / 2)^2\right\}_{i \in \mathbb{N}}\) and \(\left\{(i \pi / 2)^2\right\}_{2 i \in \mathbb{N}}\) with corresponding eigenfunctions \(\phi_0\), \(\left\{\phi_i\right\}_{i \in \mathbb{N}}\), and \(\left\{\psi_i\right\}_{2 i \in \mathbb{N}}\) given by \(\phi_0(s)=1 / \sqrt{3}\) and \[\begin{equation*} \phi_i(s)=C_{\phi, i}\begin{cases} -2 \sin (\frac{i\pi}{2}) \cos (\frac{i \pi t}{2}), & s \in e_1, \\ \sin (i \pi t / 2), & s \in e_2, \end{cases}, \quad \psi_i(s)=\frac{\sqrt{3}}{\sqrt{2}} \begin{cases} (-1)^{i / 2} \cos (\frac{i \pi t}{2}), & s \in e_1, \\ \cos (\frac{i \pi t}{2}), & s \in e_2, \end{cases}, \end{equation*}\] where \(C_{\phi, i}=1\) if \(i\) is even and \(C_{\phi, i}=1 / \sqrt{3}\) otherwise. Moreover, these functions form an orthonormal basis for \(L_2(\Gamma_T)\).
graph <- gets.graph.tadpole(h = 0.01)
eigenpairs <- gets.eigen.params(N_finite = 25, kappa = 1, alpha = 0.5, graph)
EIGENFUN <- eigenpairs$EIGENFUN
EIGENVAL <- eigenpairs$EIGENVAL
INDEX <- eigenpairs$INDEXFigure 2: Illustration of the eigenfunctions on the tadpole graph.
The piecewise constant projection \(U_{\text{piecewise}}(t^*_\ell)\) of the approximated values \(U_{\text{approx}}(t_k)\), defined on a coarse time grid \(\{t_k\}_{k=0}^{N}\), onto a finer grid \(\{t^*_\ell\}_{\ell=0}^{M}\), is given by
\[\begin{equation} U_{\text{piecewise}}(t^*_\ell) = \begin{cases} U_{\text{approx}}(t_0), & \text{if } t^*_\ell = 0 \\\\ U_{\text{approx}}(t_k), & \text{if } t^*_\ell \in (t_{k-1}, t_k], \quad \text{for } k = 1, \dots, N. \end{cases} \end{equation}\]
This defines a function that is constant on each interval \((t_{k-1}, t_k]\), and takes the value of the approximation at the right endpoint \(t_k\) of the interval.
See the Fuctionality
page for the implementation of the function
construct_piecewise_projection() that performs this
operation.
A function \(U_{\text{coarse}}(s)\) defined on a coarse mesh with \(N_{h}\) nodes can be projected onto a fine mesh with \(N_{h_{\text{ok}}}\) nodes by doing \(U_{\text{fine}}(s) = \boldsymbol{\Psi} U_{\text{coarse}}(s)\), where \(\boldsymbol{\Psi}\) is a matrix with entries \(\boldsymbol{\Psi}_{ij}=\psi^j_h(s_i)\) for \(j=1\dots, N_h\) and \(i = 1,\dots N_{h_{\text{ok}}}\).
# Parameters
T_final <- 2
kappa <- 15
N_finite = 4 # choose even
adjusted_N_finite <- N_finite + N_finite/2 + 1
# Coefficients for u_0 and f
coeff_U_0 <- 50*(1:adjusted_N_finite)^-1
coeff_U_0[-5] <- 0
coeff_FF <- rep(0, adjusted_N_finite)
coeff_FF[7] <- 10
AAA = 1
OMEGA = pi
# Overkill parameters
overkill_time_step <- 0.01
overkill_h <- 0.01
# Finest time and space mesh
overkill_time_seq <- seq(0, T_final, length.out = ((T_final - 0) / overkill_time_step + 1))
overkill_graph <- gets.graph.tadpole(h = overkill_h)
# Compute the weights on the finest mesh
overkill_graph$compute_fem() # This is needed to compute the weights
overkill_weights <- overkill_graph$mesh$weights
alpha <- 0.5 # from 0.5 to 2
beta <- alpha / 2
# Compute the eigenvalues and eigenfunctions on the finest mesh
overkill_eigen_params <- gets.eigen.params(N_finite = N_finite, kappa = kappa, alpha = alpha, graph = overkill_graph)
EIGENVAL_ALPHA <- overkill_eigen_params$EIGENVAL_ALPHA # Eigenvalues (they are independent of the meshes)
overkill_EIGENFUN <- overkill_eigen_params$EIGENFUN # Eigenfunctions on the finest mesh
# Compute the true solution on the finest mesh
overkill_U_true <- overkill_EIGENFUN %*%
outer(1:length(coeff_U_0),
1:length(overkill_time_seq),
function(i, j) (coeff_U_0[i] + coeff_FF[i] * G_sin(t = overkill_time_seq[j], A = AAA, lambda_j_alpha_half = EIGENVAL_ALPHA[i], omega = OMEGA)) * exp(-EIGENVAL_ALPHA[i] * overkill_time_seq[j]))
h <- 0.2
time_step <- 0.2
m <- 1
graph <- gets.graph.tadpole(h = h)
graph$compute_fem()
G <- graph$mesh$G
C <- graph$mesh$C
L <- kappa^2*C + G
eigen_params <- gets.eigen.params(N_finite = N_finite, kappa = kappa, alpha = alpha, graph = graph)
EIGENFUN <- eigen_params$EIGENFUN
U_0 <- EIGENFUN %*% coeff_U_0 # Compute U_0 on the current mesh
A <- graph$fem_basis(overkill_graph$get_mesh_locations())
time_seq <- seq(0, T_final, length.out = ((T_final - 0) / time_step + 1))
my_op_frac <- my.fractional.operators.frac(L, beta, C, scale.factor = kappa^2, m = m, time_step)
INT_BASIS_EIGEN <- t(overkill_EIGENFUN) %*% overkill_graph$mesh$C %*% A
# Compute matrix F with columns F^k
FF_approx <- t(INT_BASIS_EIGEN) %*%
outer(1:length(coeff_FF),
1:length(time_seq),
function(i, j) coeff_FF[i] * g_sin(r = time_seq[j], A = AAA, omega = OMEGA))
U_approx <- matrix(NA, nrow = nrow(C), ncol = length(time_seq))
U_approx[, 1] <- U_0
for (k in 1:(length(time_seq) - 1)) {
U_approx[, k + 1] <- as.matrix(my.solver.frac(my_op_frac, my_op_frac$C %*% U_approx[, k] + time_step * FF_approx[, k + 1]))
}
projected_U_approx <- A %*% U_approx
projected_U_piecewise <- construct_piecewise_projection(projected_U_approx, time_seq, overkill_time_seq)graph.plotter.3d(overkill_graph, overkill_time_seq, overkill_time_seq, overkill_U_true, projected_U_piecewise)